*
* $Id$
*
* $Log: gphxsi.F,v $
* Revision 1.2  2002/12/02 16:37:45  brun
* Changes from Federico Carminati and Peter Hristov who ported the system
* on the Ithanium processors.It is tested on HP, Sun, and Alpha, everything
* seems to work. The optimisation is switched off in case of gcc2.xx.yyy
*
* Revision 1.1.1.1  2002/07/24 15:56:26  rdm
* initial import into CVS
*
* Revision 1.1.1.1  2002/06/16 15:18:41  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:20  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:21:31  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.22  by  S.Giani
*-- Author :
      SUBROUTINE G3PHXSI
*.
*.    ******************************************************************
*.    *  Creates PHXS bank containing x-section constants              *
*.    *                                                                *
*.    *    ==>CALLED BY : G3PHYSI                                      *
*.    *       AUTHOR : J. CHWASTOWSKI                                  *
*.    *                                                                *
*.    ******************************************************************
*.
*      The photoelectric effect x-section bank
*       1                 - Number of elements neded to create medium = NZ
*       2      <-> NZ+1   -  Z of elements
*       NZ+1   <-> 2*NZ+1 -  0
*       2+2*NZ <-> 3*NZ+1 -  weight of x-section constants
*       3*NZ+2 <-> top    -  X-section constants
#include "geant321/gcbank.inc"
#include "geant321/gcphys.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcunit.inc"
#include "geant321/gcpmxz.inc"
#include "geant321/gcphxs.inc"
#include "geant321/gc10ev.inc"
      COMMON / FIXIT / JPHXS, JPHFN
      CHARACTER*1 CHSHEL, CHGROU
      DIMENSION EUP(MAXINT),TMP(MAXPOW,MAXINT),ESHL(24)
      PARAMETER (NXSB=73,NXSBF=70)
      Z = Q(JMA+7)
      IPHXSP = 0
* Check Z range of validity for Sandia parametrization
      IF(Z.GE.1.AND.Z.LE.MAXELZ) THEN
*
* Find number of elements neded to create current medium
c
* Is this medium a mixture ?
         NMIX = Q(JMA+11)
         IF(NMIX.GT.1) THEN
            NZ = 0
            JMIXT=LQ(JMA-5)
            DO 10 I = 1,NMIX
               ZCUR = Q(JMIXT+NMIX+I)
               IF(ZCUR.NE.INT(ZCUR)) THEN
*
* When Z is non integer you need 2 elements
                  NZ=NZ+2
               ELSE
                  NZ=NZ+1
               ENDIF
   10       CONTINUE
            CALL G3WORK(3*NZ+1)
            WS(1) = NZ
*
* Calculate weigths
*
            K = 1
            DO 20 I = 1,NMIX
               ZCUR = Q(JMIXT+NMIX+I)
               HZCUR = INT(ZCUR)
               WS(1+K) = HZCUR
               IF(ZCUR.NE.HZCUR) THEN
                  WS(1+2*NZ+K) = (HZCUR+1.-ZCUR)*Q(JMIXT+2*NMIX+I)
                  K = K+1
                  WS(1+K) = HZCUR+1
                  WS(1+2*NZ+K) = (ZCUR-HZCUR)*Q(JMIXT+2*NMIX+I)
               ELSE
                  WS(1+2*NZ+K) = Q(JMIXT+2*NMIX+I)
               ENDIF
               K = K+1
   20       CONTINUE
*
* Do Z values repeat ?
*
            K = NZ
            DO 40 I = 1,NZ-1
               Z1 = WS(1+I)
               IF(Z1.GT.0.0) THEN
                  DO 30 J = I+1,NZ
                     Z2 = WS(1+J)
                     IF(Z1.EQ.Z2) THEN
*                        WS(1+NZ+I) = WS(1+NZ+I)+WS(1+NZ+J)
                        WS(1+2*NZ+I) = WS(1+2*NZ+I)+WS(1+2*NZ+J)
                        WS(1+J) = -WS(1+J)
                        K = K-1
                     ENDIF
   30             CONTINUE
               ENDIF
   40       CONTINUE
* Now you can book a pht. eff. x-sec. constant bank and hang it as
* a first struc. link from JPHOT.
* From this bank you hang the banks for each separate Z !!!
            NW = NZ*NXSB+2
            IF(K.EQ.1) NW = NXSB+2
            CALL MZBOOK(IXCONS,JPHXS,JPHOT,-1,'PHXS',NZ,NZ,NW,3,0)
* Fill JPHXS bank and calculates the weights
            Q(JPHXS+1) = K
            NZN = K
            K = 1
            DO 50 I = 1,NZ
               IF(WS(1+I).GT.0.0) THEN
                  Q(JPHXS+1+K) = WS(1+I)
*                  Q(JPHXS+1+K+NZN) = WS(1+I+NZ)
                  Q(JPHXS+1+K+2*NZN) = WS(1+I+2*NZ)
                  K = K+1
               ENDIF
   50       CONTINUE
         ELSE
* Current medium consists of one "element"
            NZ = 1
            ZCUR = Z
            HZCUR = INT(ZCUR)
            IF(MOD(ZCUR,HZCUR).EQ.0) THEN
* Now you can book a pht. eff. x-sec. constant bank and hang it as
* a first struc. link from JPHOT.
               NW = NXSB+2
               CALL MZBOOK(IXCONS,JPHXS,JPHOT,-1,'PHXS',1,1,NW,3,0)
               Q(JPHXS+1) = 1.
               Q(JPHXS+2) = ZCUR
               Q(JPHXS+4) = 1.
            ELSE
* Somebody is cheating. We need two elements
               NZ = NZ+1
               NW = 2*NXSB+2
               CALL MZBOOK(IXCONS,JPHXS,JPHOT,-1,'PHXS',2,2,NW,3,0)
               Q(JPHXS+1) = 2.
               Q(JPHXS+2) = HZCUR
               Q(JPHXS+3) = HZCUR+1
               Q(JPHXS+6) = (HZCUR+1.-ZCUR)
               Q(JPHXS+7) = (ZCUR-HZCUR)
            ENDIF
         ENDIF
* We passed the bank booking phase and we can start real work
         NZ = Q(JPHXS+1)
* Create a temporary working space
         IBINS = 5*NZ+NZ*(MAXPOW+1)*(MAXINT+1)
         CALL G3WORK(IBINS)
         DO 60 I = 1,IBINS
            WS(I) = 0.0
   60    CONTINUE
         KS = 5*NZ
         L = 0
         DO 170 N = 1,NZ
            JPHXS = LQ(JPHOT-1)
            ZCUR = Q(JPHXS+1+N)
            IZ = ZCUR
            DO 70 I = 1,24
               ESHL(I) = 0.0
   70       CONTINUE
            CALL G3FSHLS(ZCUR,ESHL,NSHL)
            DO 80 I = 1,NSHL
               ESHL(I) = ESHL(I)*1.E-3
   80       CONTINUE
* Use Sandia data.
* Find out the interval upper limit.
            IMAX = 0
            DO 90 I = 1,MAXINT
               CHSHEL=CRNGUP(I,IZ)(1:1)
               CHGROU=CRNGUP(I,IZ)(2:2)
               IF(CHSHEL.EQ.'I') THEN
                  EUP(I) = 55.E11
                  IMAX = I
                  GO TO 100
               ELSEIF(CHSHEL.EQ.'K') THEN
                  EUP(I) = ESHL(1)
               ELSEIF(CHSHEL.EQ.'L') THEN
                  IF(CHGROU.EQ.'1') THEN
                     EUP(I) = ESHL(2)
                  ELSEIF(CHGROU.EQ.'2') THEN
                     EUP(I) = ESHL(3)
                  ELSEIF(CHGROU.EQ.'3') THEN
                     EUP(I) = ESHL(4)
                  ELSE
                     CHMAIL = ' GPHXSI Error: Inconsistent data for L '
     +               //'shell.  Maximum coded L3. Found shell '
     +               //'name: '//CRNGUP(I,IZ)
                     CALL GMAIL(0,0)
                  ENDIF
               ELSEIF(CHSHEL.EQ.'M') THEN
                  IF(CHGROU.EQ.'1') THEN
                     EUP(I) = ESHL(5)
                  ELSEIF(CHGROU.EQ.'2') THEN
                     EUP(I) = ESHL(6)
                  ELSEIF(CHGROU.EQ.'3') THEN
                     EUP(I) = ESHL(7)
                  ELSEIF(CHGROU.EQ.'4') THEN
                     EUP(I) = ESHL(8)
                  ELSEIF(CHGROU.EQ.'5') THEN
                     EUP(I) = ESHL(9)
                  ELSE
                     CHMAIL = ' GPHXSI Error: Inconsistent data for M '
     +               //'shell.  Maximum coded M5. Found shell '
     +               //'name: '//CRNGUP(I,IZ)
                     CALL GMAIL(0,0)
                  ENDIF
               ELSEIF(CHSHEL.EQ.'N') THEN
                  IF(CHGROU.EQ.'1') THEN
                     EUP(I) = ESHL(10)
                  ELSEIF(CHGROU.EQ.'2') THEN
                     EUP(I) = ESHL(11)
                  ELSEIF(CHGROU.EQ.'3') THEN
                     EUP(I) = ESHL(12)
                  ELSEIF(CHGROU.EQ.'4') THEN
                     EUP(I) = ESHL(13)
                  ELSEIF(CHGROU.EQ.'5') THEN
                     EUP(I) = ESHL(14)
                  ELSEIF(CHGROU.EQ.'6') THEN
                     EUP(I) = ESHL(15)
                  ELSEIF(CHGROU.EQ.'7') THEN
                     EUP(I) = ESHL(16)
                  ELSE
                     CHMAIL = ' GPHXSI Error: Inconsistent data for N '
     +               //'shell. Maximum coded N7. Found shell name:'
     +               //' '//CRNGUP(I,IZ)
                     CALL GMAIL(0,0)
                  ENDIF
               ELSEIF(CHSHEL.EQ.'O') THEN
                  IF(CHGROU.EQ.'1') THEN
                     EUP(I) = ESHL(17)
                  ELSEIF(CHGROU.EQ.'2') THEN
                     EUP(I) = ESHL(18)
                  ELSEIF(CHGROU.EQ.'3') THEN
                     EUP(I) = ESHL(19)
                  ELSEIF(CHGROU.EQ.'4') THEN
                     EUP(I) = ESHL(20)
                  ELSEIF(CHGROU.EQ.'5') THEN
                     EUP(I) = ESHL(21)
                  ELSE
                     CHMAIL = ' GPHXSI Error: Inconsistent data for O '
     +               //'shell.  Maximum code O5. Found shell name:'
     +               //' '//CRNGUP(I,IZ)
                     CALL GMAIL(0,0)
                  ENDIF
               ELSEIF(CHSHEL.EQ.'P') THEN
                  IF(CHGROU.EQ.'1') THEN
                     EUP(I) = ESHL(22)
                  ELSEIF(CHGROU.EQ.'2') THEN
                     EUP(I) = ESHL(23)
                  ELSEIF(CHGROU.EQ.'3') THEN
                     EUP(I) = ESHL(24)
                  ELSE
                     CHMAIL = ' GPHXSI Error: Inconsistent data for P '
     +               //'shell.  Maximum coded P3. Found shell '
     +               //'name: '//CRNGUP(I,IZ)
                     CALL GMAIL(0,0)
                  ENDIF
               ELSE
                  READ(CRNGUP(I,IZ),'(F6.0)') EUP(I)
               ENDIF
               IF(EUP(I).EQ.0.0) THEN
                  WRITE(CHMAIL,'(A44,I5)') ' GPHXSI Error: Upper limit '
     +            //'= 0. Interval:',I
                  CALL GMAIL(0,0)
                  STOP 14
               ENDIF
   90       CONTINUE
  100       CONTINUE
            DO 120 I = 1,IMAX
               DO 110 J = 1,MAXPOW
                  TMP(J,I) = COFS(J,I,IZ)*Q(JPHXS+1+2*NZ+N)
  110          CONTINUE
  120       CONTINUE
* Copy upper limits and coofficients to a work bank
            K = KS+1
            WS(K) = GPOMIN(IZ)
            WS(K+1) = 0.0
            WS(K+2) = 0.0
            WS(K+3) = 0.0
            WS(K+4) = 0.0
            K = K+5
            DO 140 I = 1,IMAX
               WS(K) = EUP(I)
               K = K+1
               DO 130 J = 1,MAXPOW
                  WS(K) = TMP(J,I)
                  K = K+1
  130          CONTINUE
  140       CONTINUE
            KS = KS+(MAXINT+1)*(MAXPOW+1)
            IWS(NZ+N) = IMAX
            IWS(2*NZ+N) = 0
            IWS(3*NZ+N) = 0
*            WS(N) = Q(JPHXS+1+2*NZ+N)
* Create element x-section & final state bank
            CALL MZBOOK(IXCONS,JPHFN,JPHXS,-N,'PHFN',0,0,IMAX*5+1,3,0)
            Q(JPHFN+1) = IMAX
* Update pointer
            JPHFN = JPHFN+1
            KFN = 1
* Copy energy & x-section parameters
            DO 160 I = 1,IMAX
               Q(JPHFN+KFN) = EUP(I)
               KFN = KFN+1
               DO 150 J = 1,MAXPOW
c                  Q(JPHFN+KFN) = TMP(J,I)*Q(JPHXS+1+2*NZ+N)
                  Q(JPHFN+KFN) = TMP(J,I)
                  KFN = KFN+1
  150          CONTINUE
  160       CONTINUE
* Get shells decay parameters
            CALL G3FSHDC(N,ZCUR)
  170    CONTINUE
* Now find intervals and calculate the coofs for each
         K = 0
         JPHXS = LQ(JPHOT-1)
         IF(NZ.LT.2) THEN
* Simple element so life is easy and nice
            Q(JPHXS+5) = IMAX+1
            JPHXS6 = JPHXS+6
            Q(JPHXS6) = GPOMIN(IZ)
            Q(JPHXS6+1) = 0.0
            Q(JPHXS6+2) = 0.0
            Q(JPHXS6+3) = 0.0
            Q(JPHXS6+4) = 0.0
            JPHXS6 = JPHXS6+5
            DO 190 I = 1,IMAX
               Q(JPHXS6+K) = EUP(I)
               K = K+1
               DO 180 J = 1,MAXPOW
                  Q(JPHXS6+K) = TMP(J,I)
                  K = K+1
  180          CONTINUE
  190       CONTINUE
            NPSH = NW-5*(IMAX+1)-5
         ELSE
* More elements. It will not be so easy
* The following code is difficult and probably there are better solutions
* but I could think only about this one.
            IPHXSP = JPHXS+2+3*NZ
            IPIMAX = JPHXS+2+3*NZ
            IOFFST = (MAXPOW+1)*(MAXINT+1)
            DO 250 II = 1,NZ*(MAXINT+1)
* Find the interval for which the upper limit is the smallest
               AMINV = 1.E20
               DO 200 I = 1,NZ
                  K = IWS(2*NZ+I)
                  IPOINT = 5*NZ+1+(I-1)*IOFFST+K
                  IEUP = 4*NZ+I
                  WS(IEUP) = WS(IPOINT)
                  IF(WS(IEUP).LT.AMINV) AMINV = WS(IEUP)
  200          CONTINUE
               L = 0
               DO 210 I = 1,NZ
                  IF(WS(4*NZ+I).LE.AMINV) L = L+1
  210          CONTINUE
               IF(L.LT.1) THEN
                  CHMAIL = ' GPHXSI Error: Zero intervals found.'
                  CALL GMAIL(0,0)
                  STOP 16
               ENDIF
* Copy to JPHXS bank
               Q(IPIMAX) = Q(IPIMAX)+1.
               Q(IPHXSP+1) = AMINV
               IPHXSP = IPHXSP+1
               DO 230 J = 1,MAXPOW
                  QS = 0.0
                  DO 220 I = 1,NZ
                     K = IWS(2*NZ+I)
                     IPOINT = 5*NZ+1+(I-1)*IOFFST+K+J
c                     QS = QS+WS(I)*WS(IPOINT)
                     QS = QS+WS(IPOINT)
  220             CONTINUE
                  Q(IPHXSP+1) = QS
                  IPHXSP = IPHXSP+1
  230          CONTINUE
               IF(L.EQ.NZ.AND.AMINV.EQ.55.E11) GO TO 260
* Update local pointers
               DO 240 I = 1,NZ
                  IEUP = 4*NZ+I
                  IF(WS(IEUP).LE.AMINV) THEN
                     INZ2 = 2*NZ+I
                     IF(IWS(3*NZ+I).LT.IWS(I+NZ)) THEN
                        IWS(INZ2) = IWS(INZ2)+5
                        IWS(3*NZ+I) = IWS(3*NZ+I)+1
                     ENDIF
                  ENDIF
  240          CONTINUE
  250       CONTINUE
* It is THE END of th x-secs. part, however you may not believe it.
  260       CONTINUE
            NIT=Q(IPIMAX)
            IIT = -1
  261       IIT = IIT+1
*
* It just may happen that some of the energy limits are the same
               IENE=IIT*5+1
               IF(ABS(Q(IPIMAX+IENE)-Q(IPIMAX+IENE+5)).LT.
     +                                Q(IPIMAX+IENE)*5E-4) THEN
                  DO 262 II=1,(NIT-IIT-1)*5
                     Q(IPIMAX+IENE+II)=Q(IPIMAX+IENE+II+5)
  262             CONTINUE
                  NIT=NIT-1
                  IIT=IIT-1
               ENDIF
            IF(IIT.LT.NIT-1) GO TO 261
            Q(IPIMAX)=NIT
            NPSH = NW-Q(IPIMAX)*5-3*NZ-2
         ENDIF
* Return unused locations in JPHXSI bank to the system
         IF(NPSH.GT.0) THEN
            JPHXS = LQ(JPHOT-1)
            CALL MZPUSH(IXCONS,JPHXS,0,-NPSH,'R')
         ENDIF
      ELSEIF(Z.GT.100.) THEN
* Just in case we got called
         CHMAIL = ' GPHXSI Error: Z > 100. No Sandia parameters. '
         CALL GMAIL(0,0)
      ENDIF
      END
